Statistical-mechanical iterative algorithms on complex networks 
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The Ising models have been appUed for various problems on information sciences, social sciences, 
and so on. In many cases, solving these problems corresponds to minimizing the Bethe free energy. 
To minimize the Bethe free energy, a statistical- mechanical iterative algorithm is often used. We 
study the statistical-mechanical iterative algorithm on complex networks. To investigate effects of 
heterogeneous structures on the iterative algorithm, we introduce an iterative algorithm based on 
information of heterogeneity of complex networks, in which higher-degree nodes are likely to be 
updated more frequently than lower-degree ones. Numerical experiments clarified that the usage 
of the information of heterogeneity affects the algorithm in Barabasi and Albert networks, but 
does not influence that in Erdos and Renyi networks. It is revealed that information of the whole 
system propagates rapidly through such high-degree nodes in the case of Barabasi- Albert's scale-free 
networks. 
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I. INTRODUCTION 



In recent years, complex networks have been studied a 
lot because of their importance as backbones of real com- 
plex systems iljj, l^j- Many models have been proposed 
to capture properties of real-world networks; the rep- 
resentatives of them are small-world networks proposed 
by Watts and Strogatz and scale- free networks pro- 
posed by Barabasi and Albert Q . Much elaboration has 
made clear their topological features such as the small- 
world property, the clustering property, and community 
structures. Not only their topological features, but also 
their dynamical properties have been clarified. For exam- 
ple, it has been revealed that the scale-free structure has 
large influence on dynamics in its heterogeneous struc- 
ture ^ Jj. It has been shown that the epidemic 
threshold vanishes if their degree distributions have a 
power-law form, P{k) ~ k~'^, where the characteristic 
degree exponent 7 is lower than 3 |^. Other dynamics 
are also important in order to understand characteris- 
tic properties of heterogeneous structures. Indeed, there 
are many examples of dynamics on complex networks, in 
which nodes are interacting each other and a macroscop- 
ically functional behavior emerges. Examples of these 
networks with dynamics are neural networks in a brain 
0, Hopfield models in complex networks |^, 0, , pro- 
tein folding ^12.] , and combinatoric optimization problems 

While the heterogeneous structures of the complex net- 
works have large influence on their dynamics, one can use 
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efficiently the heterogeneous property to control their dy- 
namical properties. The usage of the heterogeneous prop- 
erties has been studied a lot mainly from the viewpoint of 
statistical physics. For instance, in epidemic dynamics, 
targeted immunization schemes based on the degree hier- 
archy work efficiently, which recover the epidemic thresh- 
old in a case of the degree distribution P{k) ~ k~'', 
where 7 < 3 0|. Therefore, one expects that high- 
degree nodes, so-called hubs, play important roles in the 
dynamics on the heterogeneous networks, and hence, one 
could use the heterogeneity to make efficient algorithms 
in a wide variety of problems. This is true in some cases, 
and an efficient local search algorithm has been proposed, 
which efficiently uses information of high-degree nodes 
(hubs) 

The Ising spin system is one of the simple dynam- 
ical models, and has attracted a lot of attentions by 
their theoretical interests. In recent years, their theoret- 
ical analyses on complex networks have been performed 
[TM [T9I I20I I21I I . The Ising spin system could represent 
various real world phenomena; simple opinion dynamics 
in a society is described by two Ising spin states. More- 
over, the Ising spin system has been used for many prob- 
lems on information sciences such as an artificial intelli- 
gence, image processing, error-correcting codes, and opti- 
mization problems . These problems can be solved by 
minimizing the free energy with suitably selected Hamil- 
tonian. Though it is difhcult to obtain exact solutions 
of these problems, there are many cases in which ap- 
proximate solutions by means of the Bethe approxima- 
tion are sufficient to solve these problems. A statistical- 
mechanical iterative algorithm, so-called the Belief Prop- 
agation method _23, J4j 25}^ enables us to solve these 
problems. The iterative algorithm yields solutions ap- 
proximately equivalent to those by means of the Bethe 
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approximation. It has been clear that the iterative al- 
gorithm is useful to solve many problems on information 
sciences [2^ . 

Recently, Mooij and Kappen have studied the validity 
of the Bethe approximation |2^ They have con- 

cluded that the Bethe approximation is more powerful 
than the naive mean field approximation. A relationship 
between the convergence property of the iterative algo- 
rithm and the spin-glass phase transition has also been 
discussed, and it has been revealed that the iterative algo- 
rithm does not work enough well in the spin-glass phase. 
However, effects of heterogeneity of complex networks on 
the iterative algorithm have not been discussed yet. 

In the present paper, we discuss effects of heterogeneity 
of complex networks on the iterative algorithm. To inves- 
tigate the effects, wc introduce a statistical-mechanical it- 
erative algorithm in which higher-degree nodes are likely 
to be updated more frequently than lower-degree ones. 
Applying the newly proposed iterative algorithms to or- 
dinary random networks and scale-free networks, we re- 
veal that information of the whole systems propagates 
rapidly through the high-degree nodes; the effects enable 
the iterative algorithm to work efficiently. 

The outline of the present paper is as follows. In 
Sec. II, we introduce an Ising spin system. The Bethe 
free energy to be minimized is also shown in Sec. II. The 
statistical-mechanical iterative algorithms are explained 
in Sec. III. We show numerical results in Sec. IV. The nu- 
merical results show that the influence of heterogeneity 
appears remarkably in the case of the scale- free networks. 
Finally, we draw the main conclusions in Sec. V. 



II. THE ISING SPIN SYSTEM 

We explain an Ising spin model used in the present 
paper. Mooij and Kappen |2^ have discussed the 
spin-glass phase transitions using the similar model. 

We describe a graph G — {Q, E) as an undirected 
graph without closed one-edge loops and multiple edges, 
where Vt = {l,2,...,iV} is a set of nodes and E C 
{(ij j) 1 1 < * < i < N} is a set of edges. We denote a set 
of neighbors on node i as di. Because there are no closed 
one-edge loops and multiple edges in the graph G, the ad- 
jacency matrix A = {Aij \ i,j G J7} is defined as follows: 
Aij = Aji = 1 if e E, and Aij = Aji = other- 

wise. The degree of node i, ki, represents that node i is 
connected to ki other nodes; hence ki — \di\ — X^jGO ■ 

On the graph G, we consider a random Ising model 
which is described by the Hamiltonian 

H=- ^ J^JS^Sj, (1) 

(i,i)eB 

where Si represents an Ising spin on node i and takes 
the value -1-1 or —1, and Jij is an interaction between 
nodes i and j. Each interaction Jij is different from each 
other. We here denote the probability density function 



of Jij as Pj{Jij). From the Hamiltonian of Eq. (Q, we 
calculate the corresponding Boltzmann distribution over 
the configurations, s = (si, S2, . . . , sjy) G { — 1, +1}^, as 
follows: 

P(s) = le-'^^ = i n cb..As^,s,l (2) 

where Z is the partition function (or one may say that 
Z is the normalization constant), and f3 the inverse tem- 
perature. The notation <j)ij{si,Sj) = exp (f3JijSiSj) is 
introduced for simplicity. 

Next, we write the free energy of the Hamiltonian of 
Eq. ^ by means of the Bethe approximation. To obtain 
the Bethe free energy, the cluster variation method is 
used, in which the Bethe free energy is represented as 
a function of the marginal probability distributions, i.e., 
Pij{si,Sj) and Pi{si). As a result, the Bethe free energy 
is written as follows [25,28]: 

PJ^Bcthc{{Pi{Si)}, {Pij{Si, Sj)}) 

= - ^ ^ p,j {si , Sj ) In (j)^j (s, , Sj ) 

{i,])eE Si = ±l Sj=±l 
N 

+ Y{^-h) Y Pi{si)\npi{si) 

i=l Si = ±l 

J2 J2 P'ijisi,Sj)\npij{s,,Sj). 

(3) 

In many problems on information sciences, we require the 
solution minimizing the Bethe free energy of Eq. Q . We 
can easily see that the paramagnetic solution is a possible 
solution: 

(si) ~ J2 ''^pM) = 0, (4) 

s,=±l 

(siSj) ~ ^ ^ SjSjPy(si, Sj) = tanhJy. (5) 

Si=±l Sj=±l 

For small values of the interaction terms {Jij}, the Hes- 
sian of J-'Bethe_^ that point makes sure that the solution 
is minimum \2(i l27| . 

For simplicity, we here assume that the interactions 
{Jij} are independent Gaussian random variables with 
mean and variance 1. In this case, it has been shown 
that at the critical inverse temperature /3c, the Ising 
spin system causes the transition from the paramagnetic 
phase to a s pin-glas s phase when the inverse temperature 
(3 increases |26ll27| . 

III. STATISTICAL-MECHANICAL ITERATIVE 
ALGORITHM 

A. Update procedure 

The statistical-mechanical iterative algorithm, so- 
called Belief Propagation, is often used in various prob- 
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FIG. 1: One of examples of the update list, 
numbers in the update list in turn. 



We use node 



lems mmilll^Hlal. The iterative algorithm cal- 
culates the solution which minimizes the Bethe free en- 
ergy. Though the statistical-mechanical iterative al- 
gorithm is often performed by a synchronous update 
method, we use an asynchronous update method to in- 
vestigate effects of heterogeneous structures on the itera- 
tive algorithm. The asynchronous update method of the 
statistical-mechanical iterative algorithm needs an up- 
date list in which the node numbers to be updated are 
listed. Figure ^ shows one of the examples of the up- 
date list. In this case, at first we update the information 
related to node 4, and at the next step, update the in- 
formation related to node 26. If we reach the end of the 
update list, one iteration procedure is finished. 

The iterative algorithm by means of the Bethe approx- 
imation is as follows: 

Step 1. We set initial values of {mij{sj)\i, j G fl}. These 
initial values are chosen randomly with an uniform 
distribution (0 < mij{sj) < 1). The value mij(sj) 
is called a message from node i to node j in state 

Si. 



Step 2. We set fhij{sj) 



j{sj). The values {mijisj)} 



represent old states, which are needed to evaluate 
a convergence condition introduced later. 

Step 3. We set a node number in the update list as the value 
ofi. If it is the first step in one iteration procedure, 
we use the first node number in the update list, 
otherwise we use the next node number following 
the one used in the previous step. 

Step 4- For the nearest neighbor nodes of node i, i.e., j £ 
di, we update mij(sj) as follows: 



Y.a. (pijiSi, Sj) Ukedi\j ^kziSz) 



kedi\j 



mki{si 



(6) 



Step 5. Until we reach the end of the update list, we repeat 
Step 3 and Step 4- 

Step 6. When we reach the end of the update list, we say 
that one iteration procedure is finished. We here 
calculate the following convergence condition: 

XIXUI I'^y ~ "^iJ < (7) 



2\E\ 



If not convergent, we back to Step 2, otherwise cal- 
culate Pi{si) {i € O) using the following equation: 



P^is^) 



(8) 



(a) Method A 

I 1 I 2 I 3 I 4 I 5 I 



regular sequence 



112 13 M 15 mr- 



regular sequence 



randomized sequence (regular) 

I 4 |12| 3 |54|26| 2 |83| 9 |47| 1 1 15 | 4 1 12 | 21 1 58 1 45 | T 



(b) Method B 

I 1 I 2 I 3 I 4 I 5 I 6 I 7 I I n1 I 2 | 43 | 2 | 25 | 7 1 19 | 4 | | 2g1 



preferential sequence 



regular sequence 



randomized sequence (preferential) 

I 2 I 51 1 3 I 84 1 13 I 2 | 33 | 5 1 19 | 1 | 72 | 2 1 33 | 35 | 73 | 8 | Tf 



FIG. 2: Two methods for making the update list, (a) Or- 
dinary update list in which all node numbers are selected at 
random only two times, (b) Preferential update list in which 
high-degree nodes are likely to chosen as update nodes. This 
preferential update procedure uses the information of hetero- 
geneity of the complex network structures. 



Then, a local magnetization of Ising spin i, (si), is calcu- 
lated by 



(Sj) ^ ^ SiPi{Si). 



(9) 



Si=±l 



The pair probability distribution pij{si, Sj) is calculated 
by the similar procedure; the details of the iterative algo- 
rithm is denoted in Ref. 2^,2^ 28]. After the iterations, 
we get the solution {pi{si)} and {pij{si, sj)} which min- 
imizes the Bethe free energy of Eq. (jSJ . 



B. Making the update list 

Methods of making the update list used in the update 
procedure are arbitrary; an uniformly random update list 
is usually used, in which each node number is selected at 
random. However, we expect the following assumption: 
" Usage of heterogeneity of complex networks makes al- 
gorithms and dynamics more efficient.''^ As the simplest 
case, we here assume that high-degree nodes play an im- 
portant role in the scale-free networks; the update of the 
information of the high-degree nodes causes large effects 
on the next update step because the information of the 
high-degree nodes could be used a lot in the next update 
step. Therefore, to investigate how large the high-degree 
nodes influence on the iterative algorithm, we compare 
results obtained by two different update lists; one is made 
by the ordinary random procedure, and the other is gen- 
erated by using the information of heterogeneity of the 
network structures. 

These update lists are made by the following proce- 
dure. The making procedures are illustrated in Fig. |21 
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FIG. 3: Intuitive pictures of the update procedures, (a) Or- 
dinary update method (method A in Fig.|21) (b) Preferential 
update method (method B in Fig. |21) The messages repre- 
sented by bolder arrows tend to be updated more frequently 
than those with thinner arrows. 



Method A (ordinary update list) 

(Al) Two regular sequences are set. Each regular se- 
quence has N components, and each node number 
is listed only one time in each sequence. 

(A2) We randomly select a number from these two se- 
quences. The selected number is stocked into the 
update list. Once a number is selected, the number 
is removed from the sequences. 

(A3) We repeat the procedure (A2) until there are no 
number in the two regular sequences. Note that 
when the procedure is finished, there are 2N num- 
bers in the update list, and each node number is 
listed only two times. 

Different from method A, we make the preferentially 
selected update list using the information of the network 
structures in method B. In this case, we use the informa- 
tion of a degree of each node. 

Method B (preferential update list) 

(Bl) One regular sequence is set. The regular sequence 
has N components and each node number is listed 
only one time. Additionally, one preferential se- 
quence is set. In the preferential sequence, each 



listed number is selected with the probability pro- 
portional to its corresponding node's degree; hence, 
nodes with higher-degree tend to be selected more 
frequently than those with lower-degree. 

(B2) We select a number from these two sequences at 
random. The selected number is stocked into the 
update list. Once a number is selected, the number 
is removed from the sequences. 

(B3) We repeat the procedure (B2) until there are no 
number in the two regular sequences. Note that 
when the procedure is finished, there are 27V num- 
bers in the update list. In the update list, each 
number corresponding to each node is listed at least 
one time, and several numbers emerge a lot because 
of their high-degree. 

Intuitive pictures of effects of those methods are illus- 
trated in Fig. 121 in which the messages denoted by bolder 
arrows are updated more frequently than thinner ones. 
FigureOl^a) shows the update procedure using method A. 
In Fig.|3fa), all arrows have the same width, which repre- 
sent that we update all messages in the same frequency. 
On the other hand, in Fig. Ofb) each arrow has a dif- 
ferent thickness depending on its degree. For instance, 
node i has a lot of connections to other nodes. In the up- 
date procedure of messages related to node z, messages of 
many other nodes are needed, see Eq. jnj. Furthermore, 
the updated messages on the node i are used for the other 
update procedures on the many other nodes. Since the 
messages on high-degree node are frequently updated in 
method B, it is seems that information of local update is 
easily propagated to the whole system. Then, we expect 
that method B works more efficiently than method A for 
the heterogeneous networks like the scale-free networks. 

IV. NUMERICAL RESULTS 
A. The number of iteration steps 

To investigate effects of heterogeneous structures on 
the iterative algorithm, we performed numerical experi- 
ments. In the numerical experiments, two network struc- 
tures were used. One is the random network proposed 
by Erdos and Renyi (ER networks) , and the other is 
the scale-free network proposed by Barabasi and Albert 
(BA networks) Q . The ER networks have homogeneous 
network structures, in which each node has the similar 
degree though there are some fluctuations; the degree 
distribution of the ER networks is the Poisson distribu- 
tion. On the other hand, the BA network is the rep- 
resentative of the complex networks with heterogeneity, 
in which there are a few high-degree nodes, which is of- 
ten called hubs, and a lot of low-degree ones; the degree 
distribution is P{k) ^ k^^ . 

In all the numerical experiments, we use the networks 
with the network size N = 400, the average degree 
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FIG. 4: Results of the numerical experiments in the case of 
the ER networks. (a)The number of iteration steps needed 
for convergence as a function of the inverse temperature /3. 
(b)The percentage of convergent realizations in 200 different 
realizations. (c)The Edward- Anderson order parameter gEA. 



FIG. 5: Numerical results for the case of the BA networks. 

(a) The number of iteration steps needed for convergence. 

(b) The percentage of convergent realizations in 200 different 
realizations. (c)The Edward- Anderson order parameter qea- 



(k) = 6.0. As for the convergence condition of Eq. ||7J), 
we here set e = 10~^. If 1000 iterative procedures do not 
make the system convergent, we consider the system is 
not convergent and finish the update procedure. In addi- 
tion, as far as we checked, messages {niij} are oscillated 
in the nonconvergent realizations. We have performed 
the calculations for 200 different realizations in which 
the network structures and interaction terms {Jij} are 
different from each other. All quantities discussed below 
are averaged over only convergent realizations. 

Results for the case of the ER networks are shown in 
Fig. 21 and Fig. |S1 shows those of the BA networks. Fig- 
ures EJa) and EJa) show the number of iteration proce- 



dures needed for satisfying the convergence condition of 
Eq. 13) . The percentage of convergent realizations in 200 
different ones is shown in Figs.0{b) andlS^b). As the in- 
verse temperature /3 increases, the number of convergent 
realizations decreases at some critical inverse tempera- 
tures (3c in both cases of ER and BA networks, though 
ER and BA networks have the different critical inverse 
temperatures. We can consider the phenomenon as the 
phase transition from the paramagnetic phase to the spin- 
glass phase, as discussed by Mooij and Kappen |26l l27j . 
To confirm this, we calculate the Edwards- Anderson or- 
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der parameter 



^EA 



(10) 



iefi 



which becomes positive after the transition to the spin- 
glass phase. Figures 0Jc) and |S{c) show the Edwards- 
Anderson order parameter gsA as a function of the in- 
verse temperature /3. At the point from which the num- 
ber of convergent reahzations decreases, the Edwards- 
Anderson order qea becomes positive |31|; this suggests 
that the system is in the spin-glass state above the crit- 
ical inverse temperature /3. Note that there is not an 
external field and hence the positive value of qea does 
not indicate the ferromagnetic phase. In the spin-glass 
state, it has been shown that the statistical-mechanical 
iterative algorithm does not work enough well 26, 
While improvement of the worse convergence property of 
the statistical-mechanical iterative algorithm would be 
needed for future works, it is beyond the scope of the 
present paper. Near the critical inverse temperature /3c, 
we see remarkable improvement of the number of itera- 
tion steps needed for convergence in the case of the BA 
networks, see Fig.|SJa). Near the critical inverse temper- 
ature Pc, the number of iterative procedures needed for 
convergence of method B is 0.6 or 0.7 times as many as 
that of method A. Furthermore, the percentage of con- 
vergent realizations in the case of method B increases in 
comparison with that of method A near /3c- On the con- 
trary, for the ER random networks, both methods A and 
B give similar results because of the homogeneity of the 
ER networks, as shown in Figs.Efa) and^fb). 

Finally, we note the correctness of the solution. In 
the paramagnetic phase, all the numerical solutions yield 
(si) ~ 0. Hence, both algorithms give the correct solu- 
tion. However, it is hard to evaluate the correctness of the 
numerical solutions in the spin-glass phase; it is difficult 
to obtain the exact solution in the spin-glass phase. We 
have compared the Bethe free energies of the numerical 
solutions because a solution with a lower value of the free 
energy corresponds to more correct one, and confirmed 
that the free energies calculated by both algorithms have 
similar behavior to each other. Therefore, we consider 
that both algorithms give similar solutions while more 
detailed discussion would be needed in future research. 



B. Effect on messages 



To investigate how the proposed update algorithm ef- 
fects on the update procedure, we next investigate the 
property of message mij{sj) in Eq. ©. We here define 
the following quantities to characterize the influence on 
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FIG. 6: The averaged changes of messages over nodes with 
degree fc. (a) in the case of ER networks, (b) in the case of BA 
networks. In the both cases, we set the inverse temperature 
13 = 0.2. 



the update procedure: 

(Amy(sj))fe 



-1 



(11) 



(12) 

where Sk^^k is the Kronecker delta. The quantity Arriij 
shows the change of each message between the old it- 
erative procedure and new iterative one. The quantity 
{Amij{sj))k shows that the averaged change of messages 
over nodes with degree k. Figure |H| shows the changes 
of the messages related to the nodes with degree fc, i.e., 
{Arriij {sj))k- In these numerical experiments, we used 
the inverse temperature /3 = 0.2 in which the both sys- 
tems on the ER and the BA networks converged in all 
realizations. In the case of the ER networks, methods 
A and B have similar results as for {Amij{sj))k, though 
method A makes the change of messages a little larger, 
as shown in Fig. IHJa) . On the contrary, in the case of 
the BA networks, method B makes the changes of mes- 
sages {Amij{sj))k larger than method A for low degree 
region fc < 60. This suggests that the frequent updates 
related to high-degree nodes make the change of the mes- 
sages on low-degree nodes larger, and hence, the number 
of iteration steps decreases because the larger changes of 
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FIG. 7: Distributions of Amij{sj) at /? = 0.2. (a) and (b) 
show the cases of the ER networks and the BA networks, 
respectively. Each inset in (a) and (b) is an enlarged figure 
at large Amij{sj) region. The dot-dashed lines are guides for 
the eyes and represent P{Amij{sj)) ~ Amij{sj)~^ . 



messages cause rapid convergence to the solution. 

Distributions of the change of messages Amy(sj) are 
shown in Fig. [7Ja) and [TJb). Each inset in Fig. [7| is 
an enlarged figure at large Amij{sj) region. In the ER 
random networks, method B makes no difference from the 
case of method A. In the BA networks, as suggested in 
Fig.inib), method B makes the changes of messages larger 
than those of method A. The inset of Fig.[7|^b) shows that 
method B makes the distribution of Amy(sj) a little up 
compared to the one of method A, which suggests that 
method B work well in the case of the BA networks. It 
is valuable to note that in both cases the distribution 
of Amij{sj) shows the power-law form with the degree 
exponent —1, as shown in Figs.lTfa) and (b), while it is 
beyond the present paper to clarify why the distribution 
obeys a power-law form, so-called Zipf 's law. 

V. CONCLUSION 

In the present paper, we have investigated effects 
of heterogeneity of complex networks on statistical- 



mechanical iterative algorithms. Numerical experiments 
have shown that the usage of information of heterogene- 
ity affects the algorithm in BA networks, but does not 
influence that in ER networks. It has been revealed that 
high-degree nodes play important roles in the iterative al- 
gorithms; the usage of information of degrees makes the 
statistical-mechanical iterative algorithm more efficient 
in the case of the BA networks. At the critical inverse 
temperature, the system has transition from the para- 
magnetic phase to the spin-glass phase. Near the criti- 
cal inverse temperature f3c, the number of iteration steps 
needed for convergence decrease, and the percentage of 
convergence increases in the case of the BA networks. 
In the case of the ER networks, however, the usage of 
information of degrees does not make differences. 

The changes of messages used in the iterative algo- 
rithm have also been investigated, and it becomes clear 
that the newly proposed iterative algorithm causes larger 
changes of messages in the BA networks. This result 
agrees with intuitive pictures that the system converges 
fast when there are the larger changes of the messages. 
These results show that the information of the whole sys- 
tem propagates rapidly through high-degree nodes, and 
hence, the usage of the heterogeneity makes the iterative 
algorithm more efficient. 

We conclude that the heterogeneity is important for 
the iterative algorithms and the usage of the heterogene- 
ity makes it possible to perform the iterative calculation 
efficiently. While the current study is meaningful as a 
first step of improvement of the statistical-mechanical it- 
erative algorithm, one can consider other approaches to 
improve the iterative algorithm. We used the hetero- 
geneity of degrees of nodes in the present paper, and 
other heterogeneous property such as the load or the be- 
tweenness centrality |32j |. community structure [3^ , and 
network motifs '34'| may be available for improvement 
of the iterative algorithm. We think that the improve- 
ment of the iterative algorithm is important for large 
complex networks. In order to achieve the improvement, 
it could be helpful to research the usage of the hetero- 
geneity of network structures. Additionally, these studies 
would also be valuable to understand dynamical proper- 
ties of complex networks. 
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